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1  Objectives 

Experimental  results  performed  on  scaled  robo-insects  give  us  an  idea  of 
the  resulting  flow- field  [T|.  The  evolution  of  the  vortical  flow-held  and  the 
resulting  aerodynamic  load  time  histories  are  functions  of  various  geomet¬ 
ric  and  kinematic  parameters  of  the  system.  An  extensive  investigation  of 
the  parametric  space  using  a  grid  based  Eulerian  solver  requires  enormous 
time  and  storage,  as  the  method  needs  intricate  grids  to  resolve  the  vortex 
structures  and  the  boundary  layers.  They  also  suffer  from  numerical  diffu¬ 
sion.  To  bypass  these  issues,  we  propose  to  use  a  Lagrangian  how  solver:  it 
solves  the  incompressible  Navier-Stokes  equation  in  terms  of  vorticity.  The 
computation  is  limited  to  only  non-zero  vorticity  region,  which  increases  com¬ 
putational  efficiency,  ft  is  also  free  from  grid  based  numerical  diffusion  and 
adapts  to  complex  geometries  well. 

The  broad  objective  in  the  hrst  part  of  the  research  undertaken  is  to  use  a 
Lagrangian  particle  based  tool  and  validate  its  performance  using  earlier  ex¬ 
perimental  results  for  fundamental  happing  kinematics  like  pitch  and  plunge. 
The  specihc  aim  was  to  find  a  systematic  relationship  between  these  two  fun¬ 
damental  kinematics  at  different  parametric  regimes.  Experimental  results 
for  plunge-equivalent  pitching  cases  have  been  discussed  in  |2|.  In  this  ref¬ 
erence  comparison  between  experimental  data  and  grid-based  CFD  results 
show  poor  agreement  for  pure  pitching  cases.  The  Lagrangian  solver  would 
be  used  to  simulate  these  cases  and  the  results  will  be  compared  with  that 
presented  in  [2].  Results  for  comparison  would  be  in  terms  of  vorticity  how- 
helds  stream-wise  velocity  prohles.  Another  discrepancy  that  the  study  [2J 
reports  is,  in  the  experiments  the  heave-equivalent  pure  pitching  shows  vor¬ 
tex  street  deflection  which  was  not  present  in  its  equivalent  heaving  case. 
Thus  the  study  addresses  a  fundamental  question,  if  a  heave-equivalent  pitch 
case  [3]  in  indeed  equivalent  to  its  heave  counterpart,  or  nonlinear  effects  are 
significantly  high  not  to  make  it  so.  Based  on  our  Lagrangian  solver  results, 
these  issue  would  also  be  addressed.  This  study  would  also  consider  different 
deflection  modes  as  a  function  of  the  starting  conditions. 

In  the  second  part  of  the  research  undertaken,  the  heave  and  pitch  motion 
would  be  combined  together.  Many  invertebrates’  hovering  flight  relics  on 
this  mechanism  as  a  translation  in  the  stroke  plane  along  with  wing  rota¬ 
tion.  The  kinematics  presents  the  following  important  parameters:  reduced 
frequency,  phase  relation,  stroke  amplitude,  angle  of  rotation.  Varying  these 
parameters  gives  us  the  important  insight  on  the  flow-held  evolution  and 
aerodynamic  load  generation.  This  work  will  look  at  understanding  the  un¬ 
steady  load  generating  mechanism  for  hover.  One  of  the  main  mechanism 
used  by  insects  to  generate  large  aerodynamic  loads  is  delayed  stall.  For  a 
wing  translating  at  a  large  angle  of  attack,  a  leading  edge  vortex  is  formed  at 
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each  half-stroke  which  remains  attached  to  the  body  till  the  end  of  the  half¬ 
stroke  generating  large  suction  and  lift.  At  the  end  of  a  half-stroke,  the  wing 
rotates  and  moves  into  the  vortical  region  created  in  its  previous  half-stroke, 
sometimes  referred  to  as  ’wake  capture’  [1].  The  mechanism  of  wing  accel¬ 
eration  and  wake  capture  were  instrumental  in  increasing  the  lift  force  |4] 
in  previous  studies  on  model  wing  [5j  and  3D  computations.  However,  the 
3D  mechanism  is  different  as  there  are  spanwise  flows  to  stabilize  the  vortex, 
though  the  whole  behavior  is  Re  dependent  [6].  The  2D  behavior  is  expected 
to  be  different.  For  MAV  wing  which  operate  at  different  Reynolds  number 
ranges  than  the  natural  species,  there  is  a  need  to  revisit  the  exercises  to 
resolve  its  unsteady  flow-held.  The  Re  range  of  our  concern  is  around  1000 
and  the  how  lacks  the  stabilizing  effect  of  span-wise  how. 

It  is  also  of  interest  to  modify  the  present  code  and  make  it  more  efficient. 
The  present  code  uses  a  hrst  order  time  stepping  using  Euler’s  method.  We 
would  like  to  increase  the  time  stepping  accuracy  by  including  higher  order 
time  integration  methods  for  the  convection  step.  In  order  to  improve  the 
how-held  visualization  without  sacrificing  the  accuracy,  a  vortex  merging 
algorithm  and  some  alternate  diffusion  schemes  will  also  be  explored. 


2  Status  of  effort 

In  the  hrst  part,  the  Lagrangian  tool  has  successfully  captured  the  exper¬ 
imental  results  for  fundamental  happing  kinematics  [Hi-  The  role  of  the 
starting  condition  in  deciding  the  wake  deflection  mode  for  pure  pitch  and 
pure  plunge  is  resolved.  The  role  of  mean  angle  of  attack  has  also  been  com¬ 
mented  upon.  The  Lagrangian  tool  has  also  confirmed  the  earlier  findings 
about  the  kinematic  equivalence  of  pitch  and  plunge  for  a  symmetric  airfoil. 
It  shows  that  the  quasi-steady  criterion  and  the  Theodorsen’s  approach  give 
the  best  kinematic  equivalence  for  pitch  with  its  plunge  counterpart.  This  has 
been  investigated  for  different  stroke  amplitudes.  Simulation  has  also  been 
done  for  long  time  after  the  hrst  transients  are  settled  in  order  to  search  for 
a  second  time  stationarity.  The  observation  for  the  symmetric  prohle  is  that 
the  deflection  mode  switching  to  an  upward  pattern  irrespective  of  the  start¬ 
ing  condition  which  is  quite  similar  to  that  of  a  cambered  prohle  reported 
earlier  |8] .  We  have  published  these  results  in  \zm-  The  Lagrangian  tool  has 
also  been  modified  as  proposed  in  the  objective.  A  second  order  Runge-Kutta 
time  stepping  algorithm  has  been  successfully  implemented.  A  vortex  merg¬ 
ing  and  annihilation  scheme  has  also  been  introduced.  Both  of  them  together 
make  the  modified  tool  faster.  A  calculation  of  the  aerodynamic  loads  verify 
the  accuracy  of  the  modified  code  as  well.  A  test  case  of  a  circular  cylinder  is 
validated  successfully  with  this.  A  sinusoidal  hover  case  is  investigated  with 
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the  modified  code  at  a  moderate  Reynolds  number  around  3000.  The  flow 
field  is  compared  for  its  aerodynamic  loads  with  its  low  Reynolds  number 
counterpart  for  which  experimental  results  are  available  in  the  literature. 


3  Abstract 

A  Lagrangian  viscous  vortex  technique  is  used  in  the  present  study  to  sim¬ 
ulate  the  unsteady  flow  field  of  flapping  flight.  The  method  is  grid  free 
and  computations  are  done  only  at  the  regions  of  non-zero  vorticity.  As 
a  result  this  method  is  quite  suitable  for  simulating  the  unsteady  vortical 
flow-field  past  flapping  wings  and  airfoils;  it  provides  an  accurate  and  fast 
flow- visualization.  The  main  strength  of  this  tool  is  its  grid  free  nature;  as 
the  resolution  of  the  grid  is  a  crucial  parameter  in  resolving  the  unsteady 
flow-field  accurately.  It  is  also  relatively  easy  to  implement  a  grid-free  parti¬ 
cle  based  technique  for  a  moving  boundaries  or  complex  geometries.  In  this 
algorithm,  field  property  vorticity  is  carried  by  discrete  particles  which  con- 
vects  with  the  Biot-Savart  velocity  and  diffuses  with  random  walk.  Heave 
and  pitch  kinematics  are  fundamental  to  any  flapping  type  MAVs.  A  set 
of  pure  sinusoidal  plunging  and  its  kinematically  equivalent  pitching  using 
the  same  reduced  frequency  are  investigated.  The  focus  lies  in  investigat¬ 
ing  the  advantage  of  the  Lagrangian  method  over  grid  based  CFD  solvers 
in  capturing  the  flow  field.  Experimental  observations  have  been  taken  as 
benchmark  for  different  cases  considered.  Effect  of  mean  angle,  starting  con¬ 
ditions,  nondimensional  stroke  amplitudes  are  considered.  One  of  the  main 
questions  considered  here  is  if  there  exists  a  kinematic  equivalence  between 
sinusoidal  pitch  and  plunge. 

In  the  subsequent  part,  sinusoidal  hovering  of  a  symmetric  airfoil  at  moder¬ 
ate  Reynolds  number  is  studied.  In  the  present  study,  simulation  of  flow  over 
a  flapping  airfoil  during  sinusoidal  hover  is  studied.  Calculating  the  aerody¬ 
namic  loads  is  our  primary  concern.  The  flow  field  is  simulated  with  a  modi¬ 
fied  version  of  the  earlier  used  Lagrangian  tool.  A  higher  order  tome  stepping 
approach  is  implemented,  also  a  vortex  annihilation  and  merging  scheme  is 
used  now. The  unsteady  flow- field  is  resolved  for  a  moderate  Reynolds  num¬ 
ber  around  3000  which  is  applicable  to  various  man-made  flapping  devices 
like  Micro  Aerial  Vehicles  (MAVs).  This  flapping  behavior  is  also  compared 
with  its  low  Reynolds  number  counterpart  in  order  to  verify  any  changes  in 
aerodynamic  loads.  This  is  done  with  a  grid-based  solver  as  the  Lagrangian 
solver  shows  somewhat  poor  convergence  behavior  at  low  Reynolds  number 
regime. 
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6.2  (b)  Knowledge  resulting  form  research  in  technol¬ 

ogy  applications 

Not  applicable. 

7  Inventions 

7.1  Discoveries,  inventions,  patents 

None 

7.2  Report  of  invention  and  subcontracts 

DD  Form  882  attached. 


8  Honors/ Awards 


None 


9  Archival  Documentation 

For  the  first  part  of  the  research  effort  we  would  like  to  refer  to  our  manuscript 
accepted  for  publication  in  AIAA  journal.  We  are  attaching  a  copy  of  the 
accepted  version  of  the  manuscript.  A  technical  description  of  the  latter  part 
of  the  study  is  given  below: 

9.1  Simulation  of  Flow  over  a  Flapping  Airfoil  using  a 
Random  Vortex  Method 

9.1.1  The  modified  Lagrangian  solver 

This  is  a  grid-free  technique  and  the  flow-field  is  described  by  a  number 
of  discrete  particles  carrying  field  property  vorticity.  The  incompressible 
Navier-Stokes  equation  is  solved  in  terms  of  the  vorticity.  Velocity  field  is 
computed  from  the  Bio-Savart  law,  obtained  from  the  vector  Poisson  equation 
of  vorticity-stream  function.  Various  diffusion  models  can  be  used  in  order  to 
simulate  viscous  diffusion.  In  our  work,  a  random  walk  model  is  used  [10].  In 
some  of  our  earlier  works,  this  Lagrangian  method  was  used  to  simulate  the 
unsteady  flow-field  related  to  dynamic  stall  and  heaving-propulsion  mum 
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The  two-dimensional  incompressible  Navier-Stokes  equation  in  the  vorticity- 
transport  form  is  shown  as  [T3]: 


DQ 

~Dt 


(i) 


where  v  is  the  kinematic  viscosity,  vorticity  O  =  V  x  V  with  V  the  veloc¬ 
ity  field.  The  vorticity  and  the  velocity  field  is  given  by  a  vector  Poisson 
equation  usual 

V2!/  =  -Vx(l.  (2) 

The  solution  to  this  vector  Poisson  equation  uniquely  defines  the  velocity- 
vorticity  relationship,  known  as  the  Biot-Savart  law  Cana  IT5J  as  shown 
below. 


Pr-,0  = 


"R 


x  (r0  —  r ) 
| r0  ~  r\2 


dR  +  2 


|  r0 


■  x\ 


-dS  + 


where  fif,  is  the  rigid  body  angular  velocity  of  the  solid  whose  boundary  is 
denoted  by  S,  and  Px,  its  translational  velocity.  r0  is  the  vector  distance 
from  the  origin  of  the  reference  frame  to  the  vortex  particles  in  the  fluid 
region  R,  and  f  is  the  point  in  the  flow-field  where  the  induced  velocity  due 
to  these  vortex  particles  are  to  be  determined.  Note  that  the  velocity  field 
automatically  satisfies  the  far-held  velocity  boundary  condition  of  the  how. 
The  vorticity  held  is  discretized  as  a  collection  of  vortex  particles,  each  with 
an  associated  circulation  strength.  The  vorticity  at  any  point  in  the  how 
domain  is  given  by, 

N 

to(r)  =  '52rjf6{r-rj).  (4) 

3  = 1 

We  assume  there  ate  N  vortex  particles  at  locations  fj  with  circulations 
Tj-s  respectively;  f$  is  a  function  that  dehnes  the  vorticity  distribution,  and 
therefore  the  induced  velocity  distribution  of  the  blob.  In  the  present  work, 
Chorin  blob  model  is  used  as  the  basic  vortex  particles.  The  induced  radial 
velocity  held  due  to  a  Chorin  blob  is  given  by, 


r 

27rcr’ 

when  r  <  a 

r 

27 rr  ’ 

when  r  >  a 

(5) 

The  overall  velocity  held  consists  of  a  rotational  and  an  irrotational  part.  The 
rotational  component  of  the  how  velocity  is  obtained  from  the  velocity  in¬ 
duced  by  the  vortex  blobs.  The  contribution  from  the  solid  body  is  computed 


7 


by  approximating  the  solid  area  using  square  elements  and  summing  over  the 
contributions  from  each  of  them.  The  irrotational  component  is  obtained  by 
solving  for  potential  flow  using  the  appropriate  boundary  conditions-  the 
non-penetration  of  flow  at  the  solid  surface  in  this  case.  This  is  done  by 
using  a  panel  method  type  formulation.  According  to  this,  a  vortex  sheet 
is  assumed  to  be  attached  to  the  solid  region  boundary.  The  boundary  is 
subdivided  into  line  segments  with  a  linear  vortex  sheet  strength  distribu¬ 
tion  over  each  such  segments  or  panels.  The  induced  velocity  velocity  field  in 
complex  notation  due  to  a  panel  with  endpoints  Z\  and  z2  and  sheet  strength 
distribution  varying  from  7!  to  72  is  given  by, 


V(z)  =  u  —  iv 


where,  z'  =  (z  —  zi)eld  and  6  is  the  angle  made  by  Z\  —  z2  with  the  real 
axis.  The  velocity  distribution  over  all  panels  is  calculated  such  that  the 
condition  (^rotational  +  Kotatkmai) -n  =  Vhody.n  is  satisfied  at  specific  control 
points  along  the  surface,  which  we  choose  as  the  mid-points  of  each  panel. 
An  additional  condition  concerning  the  total  circulation  over  the  body  also 
needs  to  be  satisfied  []: 


-?l(r-  Itn 


2vr  l  VA 


Z  —  A 


+  1 


+  -{t1u 

2tt  l  A 


Z  —  A 
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+  1 


/  7 (s)ds  =  2 AbQb(t  +  At)  -  Y''  Tj(t)  (7) 

Js 

The  derived  set  of  equations  are  solved  to  obtain  the  vortex  sheet  strength 
distributed  along  the  body  surface.  This,  in  general  will  not  ensure  a  ’no¬ 
slip’  at  the  body  surface.  The  ’no-slip’  condition  14uid  =  14ody  required 
for  viscous  flows  is  ensured  by  releasing  the  Chorin  vortex  blobs  from  the 
body  surface.  These  blobs  are  generated  at  the  panel  control  points  with 
a  circulation  strength  Vsupds  and  a  core  radius  of  ds/(2ir).  Here,  ds  is  the 
individual  panel  length.  As  per  this  model,  such  a  blob  will  induce  a  core 
velocity  of  14np  in  the  tangential  direction  which  will  nullify  the  existing  slip 
velocity  at  the  body  surface. 

We  use  a  operator  splitting  technique  proposed  by  Chorin  in  which 
Eq.  ([T])  can  be  split  into  convection  and  diffusion  parts  that  are  to  be  solved 
sequentially.  This  is  represented  as, 

dfl/dt  +  V.VQ  =  0,  and  dfl/dt  =  z/V2H  (8) 

This  method  is  also  known  as  ’viscous  splitting’  and  its  convergence  has  been 
discussed  by  Beale  and  Majda  pT6] . 

The  convection  part  in  Eq.  (J8])  shows  the  invariance  of  vorticity  of  the  vortex 
particles  as  they  move  with  the  fluid.  The  convection  velocity  of  the  blobs 


are  their  self-induced  velocity  as  given  by  the  Bio-Savart  law.  A  second  order 
Runge-Kutta  time  stepping  scheme  is  employed  to  convect  the  blobs: 

Xj{t  +  At/ 2)  =  xj(t )  +  Vj[t)^- 

Xj(t  +  At)  =  Xj(t  +  At/2)  +  Vj(t  +  At/2)—.  (9) 

Here,  the  convection  displacement  of  the  jth  particle  is  denoted  by  xr 
The  solution  of  the  diffusion  part  in  Eq.  (]S])  is  given  by  a  Gaussian  probability 
density  function  [13].  A  random  walk  algorithm  ma  is  used  to  approximate 
diffusion  nans].  At  the  end  of  the  convection  step,  the  vortex  blobs  are  given 
a  diffusion  displacement  of  rjx,rjy ,  whete  rjx  and  r)y  are  random  variables  gen¬ 
erated  with  Gaussian  probability  distribution  with  zero  mean  and  standard 
deviation  of  \/2vAt.  Recently  Eldredge  [[19]  and  his  co-workers  have  also 
used  a  viscous  vortex  particle  method  (VVPM)  using  a  deterministic  diffu¬ 
sion  model  to  study  flapping  kinematics  of  airfoils.  In  another  investigation, 
A  blob  annihilation  and  merging  algorithm  is  also  used  where  blobs  of  op¬ 
posite  circulation  that  are  closer  than  a  specified  distance  are  replaced  with 
a  blob  of  circulation  strength  equal  to  the  sum  of  those  of  the  two  blobs, 
located  at  their  mid-point  and  having  their  average  blob  radius.  This  re¬ 
duces  the  computation  time  considerably  and  cleans  up  the  flow-held.  The 
aerodynamic  load  comparisons  showed  that  it  did  not  introduce  significant 
errors. 

9.1.2  Validation:  an  impulsively  started  cylinder 

The  case  of  an  impulsively  started  uniformly  moving  cylinder  has  been  stud¬ 
ied  extensively  in  the  literature  and  quite  a  number  of  excellent  benchmark 
data  are  available.  In  particular,  Koumousatkos  and  Leonard  [20]  have  ob¬ 
tained  accurate  drag  coefficient  results  from  high-resolution  computations. 
This  case  is  revisited  here  to  validate  the  present  code.  The  solid  surface  of 
the  circular  cylinder  is  subdivided  into  400  panels,  and  a  non-dimensional 
time  step  of  0.03  has  been  chosen.  Time  is  non-dimensionalized  with  respect 
to  the  cylinder  velocity  and  the  radius.  The  annihilation  distance  is  chosen 
to  be  the  average  panel  length.  The  results  obtained  for  a  Reynolds  number 
of  3000.  The  initial  flow-field  and  the  onset  of  the  Karman  vortex  structure 
is  shown  in  Fig.  [T]  and  compared  with  [21].  The  calculated  drag  coefficient 
is  seen  to  follow  the  general  trend  of  the  validation  data  [20].  Due  to  the 
stochastic  nature  of  the  random  vortex  method,  too  much  statistical  noise 
is  introduced  in  the  high  frequency  solutions.  However,  the  method  gives 
good  engineering  results  and  can  be  used  to  determine  aerodynamic  loads 
for  lower  resolution  cases.  A  slightly  lower  resolution  computation  is  also 
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carried  out  for  the  same  Reynolds  number  with  the  circular  cylinder  divided 
into  300  panels  and  a  non-dimensional  time  step  of  0.24.  The  computations 
are  carried  out  till  the  von  Karrnan  flow  patterns  is  visible  after  the  initial 
transients.  The  average  drag  coefficient  is  computed  to  be  1.49.  The  Stronhal 
number  based  on  the  frequency  of  the  vortex  shedding  is  0.2,  quite  close  to 
earlier  predicted  results  of  0.21  in  the  literature  |22j. 


Figure  1:  Starting  flow-field  for  impulsively  started  circular  cylinder  and 
validation  with  earlier  work;  left  column  is  the  present  work,  right  column  is 
from  the  earlier  computations  of  [21] 


9.1.3  Combined  Pitch  and  Plunge  in  Hover 

In  the  next  part  of  our  studies,  the  Lagrangian  code  is  modified.  A  blob 
annihilation  and  merging  algorithm  has  been  used  now.  This  makes  the 
Lagrangian  code  more  efficient  and  improves  the  flow  visualization.  Here, 
blobs  of  opposite  circulation  that  are  closer  than  a  specified  distance  are 
replaced  with  a  blob  of  circulation  strength  equal  to  the  sum  of  those  of  the 
two  blobs,  located  at  their  mid-point  and  with  their  average  blob  radius. 
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Figure  2:  Impulsively  started  circular  cylinder:  Validation  of  the  drag  coef¬ 
ficient  with  W;  solid  line  represents  the  present  work,  diamonds  are  for  the 
results  given  in  [201 


Figure  3:  Impulsively  started  circular  cylinder:  Load  time  histories. 
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This  reduces  the  computation  time  significantly  and  cleans  up  the  flow-field; 
load  comparisons  showed  that  it  did  not  introduce  significant  errors.  The 
forces  on  the  solid  body  are  now  calculated  from  the  impulse  formula  |13j: 

F  d  f  ,  dVh  A  d  ,  .  . 

—  =  — 77  /  (y  x  tydy  +  Ab— -  -  2Ab— (xb  x  ttb)  (10) 

p  at  J  D  at  at 

The  present  section  looks  at  understanding  the  unsteady  load  generation 
mechanism  during  hover.  One  of  the  main  load  generation  mechanism  for 
hovering  insect  is  delayed  stall.  For  a  wing  translating  at  a  large  angle  of 
attack,  a  leading  edge  vortex  is  formed  at  each  half-stroke  which  remains 
closely  attached  to  the  body  till  the  end  of  the  half-stroke,  generating  large 
suction  and  lift.  At  the  end  of  the  half-stroke,  the  wing  rotates  and  moves 
into  the  vortical  region  created  in  its  previous  half-stroke  and  experiences 
increased  lift,  sometimes  referred  to  as  wake  capture.  The  mechanism  of 
wing  rotation  and  wake  capture  were  found  to  be  important  contributors  to 
lift  in  previous  studies  on  model  wings  0  [23]  and  3-D  computations.  In 
3-D  flapping  wings,  spanwise  flows  stabilize  the  leading  edge  vortex,  though 
the  overall  behavior  is  dependent  on  Re  [6j.  The  2-D  behavior  is  expected 
to  be  different.  For  MAV  wings  which  operate  at  different  Re  ranges  than 
insect  species,  there  is  a  need  to  revisit  the  exercise  to  resolve  its  unsteady 
flow-field.  The  Re  range  of  our  concern  is  around  3000  and  the  flow  is  two 
dimensional. 

In  this  section,  we  present  results  for  sinusoidal  hover  m-  Once  again  a  sym¬ 
metric  airfoil  is  chosen.  The  wing  follows  a  sinusoidal  plunging  and  pitching 
motion.  Specifically  the  wing  sweeps  the  horizontal  plane  and  pitches  about 
its  spanwise  axis  with  a  frequency  /: 


x(t)  =  cos(2n ft) 

a(t )  =  «o  +  fd  sin(27 rft  +  0)  (11) 

where  x(t)  is  the  position  of  the  centre  of  the  airfoil,  and  a(t)  is  the  angle  of 
attack  with  respect  to  the  x-axis.  The  parameters  «o  ,  0,  /  and  0  are  fixed  at 
7t/2,  7t/4,  0.25  Hz  and  7t/4  .  Aq/c  is  4.0  and  the  Reynolds  number  is  chosen 
to  be  10007T,  non-dimensionalised  with  respect  to  the  maximum  translational 
velocity  of  the  centre  of  the  airfoil,  nfA0.  The  Reynolds  number  chosen  is 
much  higher  than  that  used  in  [21j  as  could  be  the  case  with  MAVs.  The 
convergence  and  accuracy  of  the  random  vortex  method  is  also  expected  to  be 
much  better  at  higher  Re  than  its  very  low  counterparts  as  the  magnitude  of 
error  in  a  random  vortex  method  is  of  0(Re -1/2)  [TO],  Larger  insects  or  MAVs 
may  operate  at  this  specified  range.  Fig JH shows  the  evolution  of  the  vorticity 
field  during  the  first  cycle  of  flapping  with  the  specified  kinematics.  Red 
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vortex  particles  show  clock-wise  circulation  and  blue  denotes  anti-clockwise 
circulation.  Since  a  stochastic  method  is  being  used  for  simulation,  the  flow 
obtained  from  different  runs  might  have  minor  variations.  But  the  major 
features  of  the  flow  always  appear. 

Due  to  the  high  angle  of  attack,  a  leading  edge  vortex  develops  that  causes 
an  initial  increase  in  lift  till  about  quarter-stroke  after  which  the  lift  begins 
to  drop  when  the  vortex  is  shed.  One  trailing  edge  vortex  is  shed  in  the 
beginning  of  the  stroke  and  by  half-stroke,  two  more  are  shed.  At  the  end  of 
the  half-stroke,  the  airfoil  changes  direction  and  moves  into  pair  of  counter 
rotating  vortices  shed  earlier.  The  airfoil  then  moves  through  them,  releasing 
two  vortex  dipoles  from  the  leading  and  trailing  edges  formed  due  to  pairing 
between  the  vortices  shed  in  the  downstroke  and  the  newly  shed  vortices  in 
the  upstroke.  The  lift  then  sees  further  increase  due  to  the  formation  of  a 
new  leading  edge  vortex,  and  decreases  again  when  it  is  shed.  At  the  end  of 
the  stroke,  the  airfoil  appears  to  move  into  another  pair  of  counter  rotating 
vortices,  but  this  time  the  trailing  edge  vortex  is  weaker.  The  variation  of 
lift  coefficient  is  plotted  in  Fig.  0  during  the  first  cycle.  The  decrease  in  lift 
observed  when  the  counter  rotating  wake  vortex  pair  is  fully  destroyed  by 
the  airfoil  moving  through  it  leads  us  to  believe  that  wake  capture  plays  an 
important  role  in  lift  generation  in  flapping  airfoils.  But  an  increase  in  lift 
is  observed  even  before  the  end  of  the  half-stroke  which  could  imply  that 
the  rotational  forces  generated  by  the  body  rotation  are  responsible  for  the 
increase  in  lift,  as  the  airfoil  has  not  moved  into  the  wake  yet.  Moreover,  the 
decrease  in  lift  also  coincides  with  the  point  of  minimum  angular  velocity. 
This  has  previously  been  dealt  with  Sun  and  Tang  [0],  whose  computations 
led  them  to  assert  that  rotational  forces  were  the  cause  as  opposed  to  the 
wake  capture  explanation  inferred  from  the  experiments  of  Dickinson  et  al  [T]. 
The  more  important  mechanism  for  lift  generation  is  the  formation  of  the 
leading  edge  vortex  when  the  airfoil  is  near  the  mean  position  of  plunging. 
The  higher  translational  velocity  as  well  as  the  delayed  stall  effect  produces 
greater  lift  than  the  aforementioned  causes  and  give  rise  to  the  higher  peak 
in  the  lift  coefficient  curve.  It  is  also  observed  that  more  vortices  shed  and 
interact  in  comparison  to  lower  Re  case.  In  particular,  dipoles  are  formed 
from  both  the  leading  and  trailing  edge  vortices  in  contrast  to  a  single  dipole 
formed  from  the  trailing  edge  alone  at  lower  Re  [25,  T9J  • 

After  the  initial  transients  it  was  observed  that  the  lift  generated  during  the 
upstroke  was  more  than  that  during  the  downstroke.  Load  time  histories  are 
plotted  in  Fig.  [6]  The  visualizations  during  the  third  and  fourth  cycle  (not 
shown  here)  reveal  that  the  leading  edge  vortex  is  not  formed  properly  during 
the  downstroke.  Hence  the  asymmetry  between  the  half-strokes.  Asymmetry 
in  the  lift  curve  between  the  half-strokes  has  also  been  observed  at  lower  Re 
for  similar  kinematics  [20], 
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Figure  4:  Evolution  of  the  vorticity  field  during  the  first  cycle  of  the  flapping 
airfoil  with  the  specified  kinematics  and  the  flow  parameters.  Red  vortex 
particles  denote  clockwise  circulation  and  blue  vortex  particles  denote  anti¬ 
clockwise  circulation.  Time  during  the  cycle  is  mentioned  as  a  fraction  of  the 
time  period. 
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Figure  5:  Variation  of  force  coefficients  over  one  cycle  for  a  sinusoidally 
flapping  wing  with  specified  kinematics  and  flow  parameters.  The  airfoil 
surface  was  divided  into  400  panels  and  the  non-dimensional  time  step  was 
O.OItt. 


Figure  6:  Variation  of  force  coefficients  over  four  cycles  for  a  sinusoidally 
flapping  wing  with  specified  kinematics  and  flow  parameters. 
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